Demonstration using the NetCDF4-Python library to access velocity data from the unstructured grid ocean model SELFE via OPeNDAP, specifying the desired URL, time, layer and lat/lon region of interest. The resulting plot of forecast velocity vectors over color-shaded bathymetry is useful for a variety of recreational and scientific purposes.
SELFE (Semi-implicit Eulerian-Lagrangian Finite Element) is an open-source community-supported modelling system, based on unstructured grids, designed for the effective simulation of 3D baroclinic circulation across river-to-ocean scales. Development is led by Dr. Joseph Zhang.
SELFE is a coupled numerical model that uses nested weather models, a coastal ocean circulation model, and a wave model. The ocean model is a triangular mesh model with horizontal resolution that is finer in complicated regions. It is layered (not depth-averaged) and includes the effects of tides, winds, and varying water densities caused by temperature and salinity changes.
In [1]:
from pylab import *
import matplotlib.tri as Tri
import netCDF4
import datetime as dt
In [2]:
# Adriatic Sea Tidal Model Data from Ivica Janekovic
url = 'http://home.irb.hr:40088/thredds/dodsC/SELFE/adria_m2_zeta_u_v_depth.nc'
nc = netCDF4.Dataset(url)
nc.variables.keys()
Out[2]:
In [3]:
# take a look at the "metadata" for the variable "u"
print nc.variables['U']
In [4]:
print nc.variables['depth']
In [5]:
# Desired time for snapshot
# ....right now (or some number of hours from now) ...
start = dt.datetime.utcnow() + dt.timedelta(hours=0)
# ... or specific time (UTC)
#start = dt.datetime(2013,3,2,15,0,0)
In [6]:
# Get desired time step
time_var = nc.variables['ocean_time']
itime = netCDF4.date2index(start,time_var,select='nearest')
# Get lon,lat coordinates for nodes (depth)
lat = nc.variables['y'][:]
lon = nc.variables['x'][:]
In [7]:
# Get Connectivity array
nv = nc.variables['ele'][:,:]-1
In [8]:
# Get water depth
h = nc.variables['h'][:]
In [9]:
dtime = netCDF4.num2date(time_var[itime],time_var.units)
daystr = dtime.strftime('%Y-%b-%d %H:%M')
print daystr
In [10]:
tri = Tri.Triangulation(lon,lat, triangles=nv)
In [11]:
# get current at layer [0 = bottom, -1 = surface]
ilayer = -1
u = nc.variables['U'][itime, :, ilayer]
v = nc.variables['V'][itime, :, ilayer]
In [12]:
levels=arange(-30,1,1) # depth contours to plot
ax= [12., 12.6, 45.2, 45.6] # region to plot
In [13]:
# find velocity points in bounding box
ind = argwhere((lon >= ax[0]) & (lon <= ax[1]) & (lat >= ax[2]) & (lat <= ax[3]))
In [14]:
# thin out the velocity vectors via random subsample
subsample=3
np.random.shuffle(ind)
Nvec = int(len(ind) / subsample)
idv = ind[:Nvec]
In [28]:
# tricontourf plot of water depth with vectors on top
fig1 = figure(figsize=(18,10))
ax1 = fig1.add_subplot(111,aspect=(1.0/cos(mean(lat)*pi/180.0)))
plt.tricontourf(tri, h,levels=levels,cmap=plt.cm.gist_earth);
plt.axis(ax)
ax1.patch.set_facecolor('0.5')
cbar=colorbar()
cbar.set_label('Water Depth (m)', rotation=-90)
Q = ax1.quiver(lon[idv],lat[idv],u[idv],v[idv],scale=2.0)
qk = quiverkey(Q,0.92,0.08,0.50,'0.5 m/s',labelpos='W')
title('SELFE Velocity, Layer %d, %s' % (ilayer, daystr));
In [38]:
# cumulative histogram of speed in displayed region
import statsmodels.api as sm
spd=sqrt(u[ind]*u[ind]+v[ind]*v[ind])
ecdf = sm.distributions.ECDF(spd.flatten())
figure(figsize=(12,4))
plot(ecdf.x,ecdf.y);grid()
title('Cumulative histogram of speed in displayed region')
xlabel('Current speed (m/s)')
Out[38]:
In [ ]: